#include"curve.h"
//用ls算法拟合函数，分为一次函数、二次函数、三次函数三个函数实现 
// 计算一次函数拟合的系数
void linearRegression(DataPoint *data, int n, double *a0, double *a1) {
    double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;
    int i;
	for (i = 0; i < n; i++) {
        sumX += data[i].x;
        sumY += data[i].y;
        sumXY += data[i].x * data[i].y;
        sumX2 += data[i].x * data[i].x;
    }
    *a1 = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);
    *a0 = (sumY - *a1 * sumX) / n;
}

// 计算二次函数拟合的系数
void quadraticRegression(DataPoint *data, int n, double *a0, double *a1, double *a2) {
    double sumX = 0, sumY = 0, sumX2 = 0, sumX3 = 0, sumX4 = 0, sumXY = 0, sumX2Y = 0;
    int i;
	for (i = 0; i < n; i++) {
        sumX += data[i].x;
        sumY += data[i].y;
        sumX2 += data[i].x * data[i].x;
        sumX3 += data[i].x * data[i].x * data[i].x;
        sumX4 += data[i].x * data[i].x * data[i].x * data[i].x;
        sumXY += data[i].x * data[i].y;
        sumX2Y += data[i].x * data[i].x * data[i].y;
    }
    double matrix[3][3] = {
        {n, sumX, sumX2},
        {sumX, sumX2, sumX3},
        {sumX2, sumX3, sumX4}
    };
    double result[3] = {sumY, sumXY, sumX2Y};
    // 这里简单使用高斯消元法求解线性方程组（实际应用可使用更优化的线性方程组解法库）
    // 下面是高斯消元的过程，将矩阵化为上三角阵
    int k,j;
	for (k = 0; k < 2; k++) {
        for (i = k + 1; i < 3; i++) {
            double factor = matrix[i][k] / matrix[k][k];
            for (j = k; j < 3; j++) {
                matrix[i][j] -= factor * matrix[k][j];
            }
            result[i] -= factor * result[k];
        }
    }
    // 回代求解
    *a2 = result[2] / matrix[2][2];
    *a1 = (result[1] - matrix[1][2] * (*a2)) / matrix[1][1];
    *a0 = (result[0] - matrix[0][1] * (*a1) - matrix[0][2] * (*a2)) / matrix[0][0];
}

// 计算三次函数拟合的系数
void cubicRegression(DataPoint *data, int n, double *a0, double *a1, double *a2, double *a3) {
    double sumX = 0, sumY = 0, sumX2 = 0, sumX3 = 0, sumX4 = 0, sumX5 = 0, sumX6 = 0;
    double sumXY = 0, sumX2Y = 0, sumX3Y = 0, sumX4Y = 0;
    int i;
	for (i = 0; i < n; i++) {
        sumX += data[i].x;
        sumY += data[i].y;
        sumX2 += data[i].x * data[i].x;
        sumX3 += data[i].x * data[i].x * data[i].x;
        sumX4 += data[i].x * data[i].x * data[i].x * data[i].x;
        sumX5 += data[i].x * data[i].x * data[i].x * data[i].x * data[i].x;
        sumX6 += data[i].x * data[i].x * data[i].x * data[i].x * data[i].x * data[i].x;
        sumXY += data[i].x * data[i].y;
        sumX2Y += data[i].x * data[i].x * data[i].y;
        sumX3Y += data[i].x * data[i].x * data[i].x * data[i].y;
        sumX4Y += data[i].x * data[i].x * data[i].x * data[i].x * data[i].y;
    }
    double matrix[4][4] = {
        {n, sumX, sumX2, sumX3},
        {sumX, sumX2, sumX3, sumX4},
        {sumX2, sumX3, sumX4, sumX5},
        {sumX3, sumX4, sumX5, sumX6}
    };
    double result[4] = {sumY, sumXY, sumX2Y, sumX3Y};
    // 高斯消元法求解线性方程组（同样可替换为更高效解法库）
    int k,j;
    for (k = 0; k < 3; k++) {
        for (i = k + 1; i < 4; i++) {
            double factor = matrix[i][k] / matrix[k][k];
            for (j = k; j < 4; j++) {
                matrix[i][j] -= factor * matrix[k][j];
            }
            result[i] -= factor * result[k];
        }
    }
    // 回代求解系数
    *a3 = result[3] / matrix[3][3];
    *a2 = (result[2] - matrix[2][3] * (*a3)) / matrix[2][2];
    *a1 = (result[1] - matrix[1][2] * (*a2) - matrix[1][3] * (*a3)) / matrix[1][1];
    *a0 = (result[0] - matrix[0][1] * (*a1) - matrix[0][2] * (*a2) - matrix[0][3] * (*a3)) / matrix[0][0];
}
